#ifdef __BORLANDC__
   #include<alloc.h>
   #include<conio.h>
#else
   #include<malloc.h>
#endif
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>

#include"draft2.h"
#include"memory2.h"
#include"error2.h"
#include"tree2.h"


/* extern  variables */
extern  StrucMesh2  mesh2;
extern  StrucTree2  tree2;


StrucDraft2  draft2;
double  SI=0.3,SM=0.2;


static void adddPoint( int v );
static void nearestDraft2( double x, double y );
static void addSize( int v );
static void remSize( int  fath );
static void nearestDraft22( double x, double y );

extern int    nTrglobal, *trglobal;
extern int    nVrglobal;
extern double *vrglobal;



void  initDraft2( void )
{
   int     i,j,v1,v2,v3;
//   char    name[128];
//   FILE    *file;
   double  x,y,s;
   long    p;

/* strcpy(name,mesh2.inFileName);   strcat(name,".vrt");
   file = fopen(name,"r");
   if( file == NULL )
     errorExit2(3,"can't open  file ");
   fscanf(file,"%d",&draft2.nPoint);
   fscanf(file,"%d",&j); 
   for(i=0;i<j;i++)
     fscanf(file,"%d",&v3); */
   draft2.nPoint = nVrglobal;
   p = (unsigned long)draft2.nPoint*sizeof(double);
   draft2.x = (double Huge *)farmalloc(p);
   draft2.y = (double Huge *)farmalloc(p);
   draft2.s = (double Huge *)farmalloc(p);
   if( draft2.x == NULL || draft2.y == NULL || draft2.s == NULL )
     errorExit2(3,"out  of memory  in  draft ");

   for(i=0;i<draft2.nPoint;i++){
/*    fscanf(file,"%lf%lf%lf",&x,&y,&s); */
      x = vrglobal[i*3];
      y = vrglobal[i*3+1];
      s = vrglobal[i*3+2];
      draft2.x[i] = x;
      draft2.y[i] = y;
      draft2.s[i] = s;
      if( s > SM )
        draft2.s[i] = SM;
   }
/* fclose(file); */

/* strcpy(name,mesh2.inFileName);   strcat(name,".tri");
   file = fopen(name,"r");
   if( file == NULL )
     errorExit2(3,"can't open  file ");
   fscanf(file,"%d",&draft2.nTria);
   fscanf(file,"%d",&j);
   for(i=0;i<j;i++)
     fscanf(file,"%d",&v3); */
   draft2.nTria = nTrglobal;
   p = (unsigned long)draft2.nTria*sizeof(int);
   draft2.v1 = (int Huge *)farmalloc(p);
   draft2.v2 = (int Huge *)farmalloc(p);
   draft2.v3 = (int Huge *)farmalloc(p);
   if( draft2.v1 == NULL || draft2.v2 == NULL || draft2.v3 == NULL )
     errorExit2(3,"out  of memory  in  draft ");

   for(i=0;i<draft2.nTria;i++){
/*    fscanf(file,"%d%d%d",&v1,&v2,&v3); */
      v1 = trglobal[i*3];
      v2 = trglobal[i*3+1];
      v3 = trglobal[i*3+2];
      draft2.v1[i] = v1-1;
      draft2.v2[i] = v2-1;
      draft2.v3[i] = v3-1;
   }
/* fclose(file); */

   draft2.neigbor = farmalloc( (MAX_NEIGBOR2+1)*sizeof(int Huge *) );
   if( draft2.neigbor == NULL )
     errorExit2(2,"farMemory");
   for(i=0;i<=MAX_NEIGBOR2;i++){
      draft2.neigbor[i] = (int Huge *) farmalloc(draft2.nPoint*sizeof(int));
      if( draft2.neigbor[i] == NULL )
        errorExit2(2,"farMemory");
   }
   for(i=0;i<draft2.nPoint;i++)
     draft2.neigbor[0][i] = 0;
   for(i=0;i<draft2.nTria;i++){
      v1 = draft2.v1[i];
      j = draft2.neigbor[0][v1];
      if( j >= MAX_NEIGBOR2 )
	errorExit2(3,"j >= NAX_NEIGBOR2");
      draft2.neigbor[0][v1]++;
      draft2.neigbor[j+1][v1] = i;

      v1 = draft2.v2[i];
      j = draft2.neigbor[0][v1];
      if( j >= MAX_NEIGBOR2 )
	errorExit2(3,"j >= NAX_NEIGBOR2");
      draft2.neigbor[0][v1]++;
      draft2.neigbor[j+1][v1] = i;

      v1 = draft2.v3[i];
      j = draft2.neigbor[0][v1];
      if( j >= MAX_NEIGBOR2 )
	errorExit2(3,"j >= NAX_NEIGBOR2");
      draft2.neigbor[0][v1]++;
      draft2.neigbor[j+1][v1] = i;
   }

   draft2.root = nearAlloc( S_StrucNode2 );
   draft2.root->flag = EMPTY;
   for(i=0;i<draft2.nPoint;i++){
      adddPoint(i);
   }

   p = (unsigned long)draft2.nPoint*sizeof(int);
   draft2.heap = nearAlloc(p);
   draft2.index= nearAlloc(p);
   draft2.nFace = 0;
   draft2.maxFace = draft2.nPoint;
   for(i=0;i<draft2.nPoint;i++){
      addSize(i);
   }
   while( draft2.nFace > 0 ){
      double  sn,dist;
      v1 = draft2.heap[0];
      s  = draft2.s[v1];
      remSize(draft2.index[v1]);
      for(i=0;i<draft2.neigbor[0][v1];i++){
         j = draft2.neigbor[i+1][v1];
         v2 = draft2.v1[j];
         if( v2 != v1 ){
            sn = draft2.s[v2];
            dist = distance(draft2.x[v1],draft2.y[v1],draft2.x[v2],draft2.y[v2]);
            if( sn-s > SI*dist ){
	       remSize(draft2.index[v2]);
               draft2.s[v2] = s+SI*dist;
               addSize(v2);
            }
         }
         v2 = draft2.v2[j];
         if( v2 != v1 ){
            sn = draft2.s[v2];
            dist = distance(draft2.x[v1],draft2.y[v1],draft2.x[v2],draft2.y[v2]);
            if( sn-s > SI*dist ){
	       remSize(draft2.index[v2]);
               draft2.s[v2] = s+SI*dist;
               addSize(v2);
            }
         }
         v2 = draft2.v3[j];
         if( v2 != v1 ){
            sn = draft2.s[v2];
            dist = distance(draft2.x[v1],draft2.y[v1],draft2.x[v2],draft2.y[v2]);
            if( sn-s > SI*dist ){
	       remSize(draft2.index[v2]);
               draft2.s[v2] = s+SI*dist;
               addSize(v2);
            }
         }
      }
   }/*draft2.nFace*/

   return;
}/*initDraft2*/


double draftSizeFace( double x, double y )
{
   int     i,v,v1=0,v2=0,v_=0,v1_=0,v2_=0,t;
   double  p,xx,x1,x2,yy,y1,y2,t1=0,t2=0,size,t1_,t2_;

   nearestDraft2(x,y);
   v = draft2.nearVert;
   for(i=0;i<draft2.neigbor[0][v];i++){
      t = draft2.neigbor[i+1][v];
      if( v == draft2.v1[t] ){
         v1 = draft2.v2[t];
         v2 = draft2.v3[t];
      }
      if( v == draft2.v2[t] ){
         v1 = draft2.v1[t];
         v2 = draft2.v3[t];
      }
      if( v == draft2.v3[t] ){
         v1 = draft2.v1[t];
         v2 = draft2.v2[t];
      }

      x1 = draft2.x[v1]-draft2.x[v];   x2 = draft2.x[v2]-draft2.x[v];
      y1 = draft2.y[v1]-draft2.y[v];   y2 = draft2.y[v2]-draft2.y[v];
      xx = x-draft2.x[v];
      yy = y-draft2.y[v];
   /* I forgot  to  write  equation */
      p = x1*y2-y1*x2;
      if( p == 0.0 )
	errorExit2(2,"imposible  p==0  in  testN");
      t1 = (y2*xx-x2*yy)/p;
      t2 =-(y1*xx-x1*yy)/p;
      if( (t1 > -1.e-10) && (t2 > -1.e-10) )
        break;/* t(v,v1,v2) is  the  number  of  triangle  ...*/
   }
   if( (t1 > -1.e-10) && (t2 > -1.e-10) )
     size = draft2.s[v]+(draft2.s[v1]-draft2.s[v])*t1+
                        (draft2.s[v2]-draft2.s[v])*t2;
   else{/* (x,y)  belong  to  empty  region  */
      nearestDraft22(x,y);
      v = draft2.nearVert;
      t1_ = t2_ = -10.;
      for(i=0;i<draft2.neigbor[0][v];i++){
         t = draft2.neigbor[i+1][v];
         if( v == draft2.v1[t] ){
            v1 = draft2.v2[t];
            v2 = draft2.v3[t];
         }
         if( v == draft2.v2[t] ){
            v1 = draft2.v1[t];
            v2 = draft2.v3[t];
         }
         if( v == draft2.v3[t] ){
            v1 = draft2.v1[t];
            v2 = draft2.v2[t];
         }

         x1 = draft2.x[v1]-draft2.x[v];   x2 = draft2.x[v2]-draft2.x[v];
         y1 = draft2.y[v1]-draft2.y[v];   y2 = draft2.y[v2]-draft2.y[v];
         xx = x-draft2.x[v];
         yy = y-draft2.y[v];
      /* I forgot  to  write  equation */
         p = x1*y2-y1*x2;
         if( p == 0.0 )
           errorExit2(2,"imposible  p==0  in  testN");
         t1 = (y2*xx-x2*yy)/p;
         t2 =-(y1*xx-x1*yy)/p;
         if( (t1 > -1.e-10) && (t2 > -1.e-10) )
           break;/* t(v,v1,v2) is  the  number  of  triangle  ...*/
         if( ( (t1 > 0 ? 0 : -t1) + (t2 > 0 ? 0 : -t2) ) <
             ( (t1_ > 0 ? 0 : -t1_) + (t2_ > 0 ? 0 : -t2_) ) ){
             t1_ = t1;
             t2_ = t2;
             v_ = v;
             v1_ = v1;
             v2_ = v2;
         }
      }
      if( (t1 > -1.e-10) && (t2 > -1.e-10) )
        size = draft2.s[v]+(draft2.s[v1]-draft2.s[v])*t1+
                           (draft2.s[v2]-draft2.s[v])*t2;
      else{
        size = draft2.s[v_]+(draft2.s[v1_]-draft2.s[v_])*t1_+
                            (draft2.s[v2_]-draft2.s[v_])*t2_;

      }
   }

   return  size;
}/*draftSizeFace*/


static void addSize( int v )
{
   int son,fath;
   int cv,*tface=draft2.heap;
   if( draft2.nFace >= draft2.maxFace )
     errorExit2(1,"draft2.nFace");

   tface[draft2.nFace] = v;  draft2.index[v] = draft2.nFace;
   son=++draft2.nFace;
aaa:fath=son/2;
   if( (fath>0) && (draft2.s[tface[son-1]] < draft2.s[tface[fath-1]])  ){
      cv = tface[son-1];
      tface[son-1] = tface[fath-1];
      tface[fath-1] = cv;
      draft2.index[cv] = fath-1;
      draft2.index[tface[son-1]] = son-1;
      son = fath;
      if( son > 1 )
        goto aaa;
   }

   return;
};/*addSize*/


static void remSize( int  fath )
{
   int son1,son2,i;
   int v,*tface;

   tface = draft2.heap;
   if( (fath >= draft2.nFace) || (fath<0) )
     errorExit2(2," draft2.nFace < 0   in  rem_face ");

   draft2.nFace--;
   tface[fath] = tface[draft2.nFace];
   draft2.index[tface[fath]] = fath;
aaa:
   i = 0;
   if( 2*fath+1 < (draft2.nFace-1) ){
      son1=2*fath+1;  son2=son1++;
      if((draft2.s[tface[fath]]<=draft2.s[tface[son1]])&&
         (draft2.s[tface[fath]]<=draft2.s[tface[son2]])) i=0;
      if((draft2.s[tface[fath]]>=draft2.s[tface[son1]])&&
         (draft2.s[tface[son1]]>=draft2.s[tface[son2]])) i=son2;
      if((draft2.s[tface[fath]]>=draft2.s[tface[son2]])&&
         (draft2.s[tface[son2]]>draft2.s[tface[son1]])) i=son1;
      if((draft2.s[tface[son1]]>=draft2.s[tface[fath]])&&
         (draft2.s[tface[fath]]>draft2.s[tface[son2]])) i=son2;
      if((draft2.s[tface[son2]]>=draft2.s[tface[fath]])&&
         (draft2.s[tface[fath]]>draft2.s[tface[son1]])) i=son1;
   }
   if( 2*fath+1 == (draft2.nFace-1) ){
      son1=2*fath+1;
      if((draft2.s[tface[fath]]>draft2.s[tface[son1]])) i=son1;else i=0;
   }

   if(i){
      v=tface[i];  tface[i]=tface[fath];   tface[fath]=v;
      draft2.index[v] = fath;
      draft2.index[tface[i]] = i;
      fath=i;   goto aaa;
   }

   return;
}; /*remSize*/


static void adddPoint( int v )
{
   int          flag,j,d,d1,bool;
   double       x,y,x1=-1.0,y1=-1.0;
   double       xc=0.5,yc=0.5,side=0.5;
   int          old_v;
   PStrucNode2  node;
   PStrucList2  list;

   node=draft2.root;   flag=node->flag;
   list = ( PStrucList2 ) node->faceNode;
   x = draft2.x[v];
   y = draft2.y[v];

   if( (x<0.)||(x>1.)||(y<0.)||(y>1) )
     errorExit2(3,"x_y_draft_insert");

   while( flag == TREE ){
      d = direction(x,y,xc,yc);  side *= 0.5;
      center(&xc,&yc,side,d);
      node = (*list)[d];
      flag = node->flag;  list = ( PStrucList2 ) node->faceNode;
   }
   if( flag == EMPTY ){
      node->flag=FACE;   node->faceNode = (PStrucFace2) v;
   }
   if( flag == FACE ){
      bool = 1;  old_v = (int) node->faceNode;
      x1 = draft2.x[old_v];
      y1 = draft2.y[old_v];
      if( (x1==x) && (y1==y) )
        errorExit2(3," coincident   in  insert  in  draft");

      while(bool){
	 node->flag = TREE;
	 node->faceNode = nearAlloc(S_StrucList2);;
	 list = ( PStrucList2 ) node->faceNode;
	 d = direction(x,y,xc,yc);
	 d1 = direction(x1,y1,xc,yc);
	 side *= 0.5;

         if( d != d1 ){
	    bool=0;
	    for(j=0;j<4;j++){
	       (*list)[j] = node = nearAlloc(S_StrucNode2);
               if(d == j){  node->flag = FACE; node->faceNode = (PStrucFace2)v;  }
               if(d1 == j){ node->flag = FACE; node->faceNode = (PStrucFace2)old_v;}
	       if( (d1 != j) && (d != j) )  node->flag = EMPTY;
	    };/* for list */
	 };/*( d != d1 )*/

	 if(d == d1){
	    for(j=0;j<4;j++){
               (*list)[j] = node = nearAlloc(S_StrucNode2);
	       if( d!=j )  node->flag = EMPTY;
	    };/* for list */
	    node = (*list)[d];  center(&xc,&yc,side,d);
	 };/*( d == d1 )*/
      };/* while  bool */
   };/* if (flag==FACE) */

   return;
};/*adddPoint*/



double distRect(double x,double y,double xe,double ye)
 /* NEED  tree2.x,tree2.y        */
{
   double p,p1;

   if( (x<=tree2.x)&&(tree2.x<=xe)&&(y<=tree2.y)&&(tree2.y<=ye) )
     return  0.;
   p  = distance(tree2.x,tree2.y,x,y);
   p1 = distance(tree2.x,tree2.y,xe,ye);  if(p1<p) p=p1;
   p1 = distance(tree2.x,tree2.y,x,ye);   if(p1<p) p=p1;
   p1 = distance(tree2.x,tree2.y,xe,y);   if(p1<p) p=p1;

   return  p;
}/*distRect*/


void nearestRec( PStrucNode2 node )
    /*  NEED    tree2.xc=0.5;tree2.yc=0.5;tree2.side=0.5;
                draft2.dist=10.0  ,tree2.x ,tree2.y
        RETURN  draft2.nearVert == nearest
    */
{
   int          flag,i,j,nn,n[4]={0,1,2,3};
   PStrucList2  list;
   int          v;
   double       x,y,s,pp,p[4];

   flag = node->flag;
   switch( flag ){
      case FACE:
         v = (int)node->faceNode;
         pp = distance(tree2.x,tree2.y,draft2.x[v],draft2.y[v]);
         if( pp < draft2.dist ){
            draft2.nearVert = v;
            draft2.dist = pp;
         }
      break;
      case TREE:
         list = ( PStrucList2 ) node->faceNode;

         x = tree2.xc;  y = tree2.yc;  s = tree2.side;
         p[0] = distRect(x-s,y-s,x,y);
         p[1] = distRect(x,y-s,x+s,y);
         p[2] = distRect(x-s,y,x,y+s);
         p[3] = distRect(x,y,x+s,y+s);

         for(j=0;j<3;j++)
           for(i=0;i<3-j;i++)
             if( p[i] > p[i+1] ){
                pp = p[i+1];  p[i+1] = p[i];  p[i] = pp;
                nn = n[i+1];  n[i+1] = n[i];  n[i] = nn;
             }

         tree2.side *= 0.5;  s *= 0.5;
         center(&tree2.xc,&tree2.yc,tree2.side,n[0]);
         if( p[0] < draft2.dist )
           nearestRec((*list)[n[0]]);

         for(i=1;i<4;i++){
            tree2.side = s;   tree2.xc = x;   tree2.yc = y;
            center(&tree2.xc,&tree2.yc,tree2.side,n[i]);
            if( p[i] < draft2.dist )
              nearestRec((*list)[n[i]]);
         }
      break;
   }/*switch*/

   return;
}/*nearestRec*/


static void nearestDraft2( double x, double y )
{
   tree2.xc=0.5;  tree2.yc=0.5;  tree2.side=0.5;
   tree2.x = x;  tree2.y = y;
   draft2.dist = 10.;
   draft2.nearVert = -1;

   nearestRec(draft2.root);
   if( draft2.nearVert == -1 )
     errorExit2(3,"draft2.nearVert == -1");

   return;
}/*nearestDraft2*/


void nearestRec22( PStrucNode2 node )
    /*  NEED    tree2.xc=0.5;tree2.yc=0.5;tree2.side=0.5;
                draft2.dist=10.0  ,tree2.x ,tree2.y
        RETURN  draft2.nearVert == nearest
    */
{
   int          flag,i,j,nn,n[4]={0,1,2,3};
   PStrucList2  list;
   int          v,v1=0,v2=0,t;
   double       x,y,s,pp,p[4];
   double       t1=0,t2=0,t1_,t2_,x1,x2,xx,y1,y2,yy;

   flag = node->flag;
   switch( flag ){
      case FACE:
         v = (int)node->faceNode;
         pp = distance(tree2.x,tree2.y,draft2.x[v],draft2.y[v]);
         if( pp < draft2.dist ){
            t1_ = t2_ = -10.;
            for(i=0;i<draft2.neigbor[0][v];i++){
               t = draft2.neigbor[i+1][v];
               if( v == draft2.v1[t] ){
                  v1 = draft2.v2[t];
                  v2 = draft2.v3[t];
               }
               if( v == draft2.v2[t] ){
                  v1 = draft2.v1[t];
                  v2 = draft2.v3[t];
               }
               if( v == draft2.v3[t] ){
                  v1 = draft2.v1[t];
                  v2 = draft2.v2[t];
               }

               x1 = draft2.x[v1]-draft2.x[v];   x2 = draft2.x[v2]-draft2.x[v];
               y1 = draft2.y[v1]-draft2.y[v];   y2 = draft2.y[v2]-draft2.y[v];
               xx = tree2.x-draft2.x[v];
               yy = tree2.y-draft2.y[v];
            /* I forgot  to  write  equation */
               pp = x1*y2-y1*x2;
               if( pp == 0.0 )
                 errorExit2(2,"imposible  p==0  in  testN");
               t1 = (y2*xx-x2*yy)/pp;
               t2 =-(y1*xx-x1*yy)/pp;
               if( (t1 > -1.e-10) && (t2 > -1.e-10) )
                 break;/* t(v,v1,v2) is  the  number  of  triangle  ...*/
               if( ( (t1 > 0 ? 0 : -t1) + (t2 > 0 ? 0 : -t2) ) <
                   ( (t1_ > 0 ? 0 : -t1_) + (t2_ > 0 ? 0 : -t2_) ) ){
                   t1_ = t1;
                   t2_ = t2;
               }
            }
            if( (t1 > -1.e-10) && (t2 > -1.e-10) ){
               draft2.nearVert = v;
               draft2.dist = pp;
            }
            else if( (t1_ > -0.15) && (t2_ > -0.15) ){
               draft2.nearVert = v;
               draft2.dist = pp;
            }
         }
      break;
      case TREE:
         list = ( PStrucList2 ) node->faceNode;

         x = tree2.xc;  y = tree2.yc;  s = tree2.side;
         p[0] = distRect(x-s,y-s,x,y);
         p[1] = distRect(x,y-s,x+s,y);
         p[2] = distRect(x-s,y,x,y+s);
         p[3] = distRect(x,y,x+s,y+s);

         for(j=0;j<3;j++)
           for(i=0;i<3-j;i++)
             if( p[i] > p[i+1] ){
                pp = p[i+1];  p[i+1] = p[i];  p[i] = pp;
                nn = n[i+1];  n[i+1] = n[i];  n[i] = nn;
             }

         tree2.side *= 0.5;  s *= 0.5;
         center(&tree2.xc,&tree2.yc,tree2.side,n[0]);
         if( p[0] < draft2.dist )
           nearestRec22((*list)[n[0]]);

         for(i=1;i<4;i++){
            tree2.side = s;   tree2.xc = x;   tree2.yc = y;
            center(&tree2.xc,&tree2.yc,tree2.side,n[i]);
            if( p[i] < draft2.dist )
              nearestRec22((*list)[n[i]]);
         }
      break;
   }/*switch*/

   return;
}/*nearestRec22*/


static void nearestDraft22( double x, double y )
{
   tree2.xc=0.5;  tree2.yc=0.5;  tree2.side=0.5;
   tree2.x = x;  tree2.y = y;
   draft2.dist = 10.;
   draft2.nearVert = -1;

   nearestRec22(draft2.root);
   if( draft2.nearVert == -1 )
     errorExit2(3,"draft2.nearVert == -1");

   return;
}/*nearestDraft22*/

